# -*- coding: utf-8 -*-
"""
Created on Tue Feb  2 15:43:45 2021

@author: rober
"""
import numpy as np

def newmarkCD(T,ksi,ag,dta,Npts):
    m = 1.0
    k = m*4*np.pi**2/T**2
    c = 4*m*ksi*np.pi/T
    dt = dta
    
    kd = 4*m/dt**2+2*c/dt+k
    
    a = np.zeros(int(Npts))
    v = np.zeros(int(Npts))
    d = np.zeros(int(Npts))
    
    for j in range(1,int(Npts)):
        dag = ag[j]-ag[j-1]
        dd = (-m*dag+2*m*a[j-1]+4*m*v[j-1]/dt+2*c*v[j-1])/kd
        dv = 2*dd/dt-2*v[j-1]
        da = 4*dd/dt**2-4*v[j-1]/dt-2*a[j-1]
        
        d[j] = d[j-1]+dd
        v[j] = v[j-1]+dv
        a[j] = a[j-1]+da
        
    Sd = max(abs(d))
    
    return Sd
